%% ex6.m
% set dir: src/
close all;
clear; clc;

METHOD = 'replicate';
PATH = '../assets/';
%% main
% load images
imgC = imread("../data/cameraman.tif");
imgE = imread("../data/einstein.tif");
imgM = rgb1gray(imread("../data/mandril_color.tif"),'NTSC');
imgL = rgb1gray(imread("../data/lena512color.tiff"),'NTSC');

% !main process
cFilt = zeros([size(imgC),6]);
eFilt = zeros([size(imgE),6]);
mFilt = zeros([size(imgM),6]);
lFilt = zeros([size(imgL),6]);

cFilt(:, :, 1) = imgC;
eFilt(:, :, 1) = imgE;
mFilt(:, :, 1) = imgM;
lFilt(:, :, 1) = imgL;
for sigma=1:5
    % compute Gaussian Kernal
    w = gaussKernel(sigma, 6 * sigma);

    % convolution
    cFilt(:, :, sigma + 1) = twodConv(imgC, w, METHOD);
    eFilt(:, :, sigma + 1) = twodConv(imgE, w, METHOD);
    mFilt(:, :, sigma + 1) = twodConv(imgM, w, METHOD);
    lFilt(:, :, sigma + 1) = twodConv(imgL, w, METHOD);
end

%% visualization
figure;
tiledlayout(4, 6, "TileSpacing", "compact");
% cameraman
nexttile;
imshow(cFilt(:, :, 1), [0, 255]);
axis off;
title('Original');
ylabel('cameraman.tif');
for i = 2:6
    nexttile;
    imshow(cFilt(:, :, i), [0, 255]);
    axis off;
    title(append('\sigma = ',string(i - 1)));
end

% einstein
nexttile;
imshow(eFilt(:, :, 1), [0, 255]);
axis off;
ylabel('einstein.tif');
for i = 2:6
    nexttile;
    imshow(eFilt(:, :, i), [0, 255]);
    axis off;
end

% mandril_color
nexttile;
imshow(mFilt(:, :, 1), [0, 255]);
axis off;
ylabel('mandril\_color.tif');
for i = 2:6
    nexttile;
    imshow(mFilt(:, :, i), [0, 255]);
    axis off;
end

% lena512color.tiff
nexttile;
imshow(lFilt(:, :, 1), [0, 255]);
axis off;
ylabel('lena512color.tiff');
for i = 2:6
    nexttile;
    imshow(lFilt(:, :, i), [0, 255]);
    axis off;
end
%saveas(gcf, append(PATH, 'ex6_results.png'), 'png')
cFiltBuildin = imgaussfilt(imgC, 1);
eFiltBuildin = imgaussfilt(imgE, 1);
lFiltBuildin = imgaussfilt(imgL, 1);
mFiltBuildin = imgaussfilt(imgM, 1);
%% build-in function
figure;
subplot(1, 4, 1);
imshow(cFiltBuildin, [0, 255]);
axis off;
subplot(1, 4, 2);
imshow(eFiltBuildin, [0, 255]);
axis off;
subplot(1, 4, 3);
imshow(lFiltBuildin, [0, 255]);
axis off;
subplot(1, 4, 4);
imshow(mFiltBuildin, [0, 255]);
axis off;
%saveas(gcf, append(PATH, 'ex6_buildin.png'), 'png')

%% diffs
cDiff = abs(cFilt(:, :, 2) - double(cFiltBuildin));
eDiff = abs(eFilt(:, :, 2) - double(eFiltBuildin));
lDiff = abs(lFilt(:, :, 2) - double(lFiltBuildin));
mDiff = abs(mFilt(:, :, 2) - double(mFiltBuildin));

figure;
tiledlayout(1, 4, "TileSpacing", "compact");
nexttile;
imshow(cDiff, [0, 2]);
axis off;
nexttile;
imshow(eDiff, [0, 2]);
axis off;
nexttile;
imshow(lDiff, [0, 2]);
axis off;
nexttile;
imshow(mDiff, [0, 2]);
axis off;

%% padding test
cFiltZero = twodConv(imgC, gaussKernel(1), 'zero');
eFiltZero = twodConv(imgE, gaussKernel(1), 'zero');
gaussKernel(1)
cDiffZR = abs(cFilt(:, :, 2) - cFiltZero);
eDiffZR = abs(eFilt(:, :, 2) - eFiltZero);

subplot(2, 1, 1);
plot(cDiffZR(50, :));
ylabel('cameraman');
xlim([-1, 900]);
subplot(2, 1, 2);
plot(eDiffZR(50, :));
ylabel('einstein');
xlim([-1, 900]);
%% rgb1gray
function g = rgb1gray(f, method)
% Transform RGB images into grayscale images using 2 optional methods.
% f(matrix): RGB image matrix.
% method(string): choose from 'average' and 'NTSC', default 'NTSC'.
if ~exist('method')
    method = 'NTSC';
end

if strcmp(method,'NTSC')
    g = f(:,:,1) * 0.2989 + f(:,:,2) * 0.5870 + f(:,:,3) * 0.1140;
elseif strcmp(method,'average')
    g = mean(f,3);
else
    error("ERROR IN PARAMS, PLEASE CHECK PARA")
end
end

%% gaussKernel
function w = gaussKernel(sig, m)
% Generate Gaussian kernel using given sigma and size.
% - Input
% sig(numeric): input grayscale matrix.
% m(numeric): convolution kernel.
% w(matrix): generated Gaussian Kernel

% If m is not provieded, calculate a properate m.
if nargin < 2
    m = ceil(6 * sig) + 1;
end

% If m is even, make it singular.
if mod(m, 2) == 0
    m = m + 1;
end

% If m is too small, throw a warning.
minSize = ceil(6 * sig) + 1;
if m < minSize
    warning('Provided m is too small.');
end

center = (m - 1) / 2;

% compute Gaussian function
w = zeros(m, m);
for i = 1:m
    for j = 1:m
        x = i - center - 1;
        y = j - center - 1;
        w(i, j) = exp(-(x^2 + y^2) / (2 * sig^2));
    end
end

% normalize Gaussian Kernel
w = w / sum(w(:));
end

%% twodConv
function g = twodConv(f,w,padding_option)
% Operate two dimensional convolution.
% - Input
% f(matrix): input grayscale matrix.
% w(matrix): convolution kernel.
% padding_option(string): padding method, providing 'replicate' and 'zero'
% option, default 'zero'

if nargin < 3
    padding_option = 'zero';
end

% determine size of kernel and image
[r, c] = size(f);
[rw, cw] = size(w);

% compute padding size
padH = floor((rw - 1) / 2);
padW = floor((cw - 1) / 2);

% padding boundary
fPadded = zeros(r + 2*padH, c + 2*padW);
switch padding_option
    case 'replicate'
        fPadded(padH + 1:end - padH, padW + 1:end - padW) = f;
        for i = 1:padH
            fPadded(i, padW + 1:end - padW) = f(1, :);
            fPadded(end - i + 1, padW + 1:end - padW) = f(end, :);
        end
        for i = 1:padW
            fPadded(padH + 1:end - padH, i) = f(:,1);
            fPadded(padH + 1:end - padH, end - i + 1) = f(:, end);
        end
        fPadded(1:padH, 1:padW) = f(1,1);
        fPadded(1:padH,end - padW + 1:end) = f(1, end);
        fPadded(end - padH + 1:end, 1:padW) = f(end, 1);
        fPadded(end - padH + 1:end, end - padW + 1:end) = f(end, end);

    case 'zero'
        fPadded(padH + 1:end - padH, padW + 1:end - padW) = f;
    otherwise
        error('Unsupported option, please check your input.')
end

g = zeros(r, c);
% conduct 2-D convolution.
for i = 1:r
    for j = 1:c
        region = fPadded(i:i + rw - 1, j:j + rw - 1);
        g(i,j) = sum(sum(region .* w));
    end
end
end
